!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==============================================================================================================================================|
subroutine weno (ro, mx, my, mz, bx, by, bz, en, nx, ng, flux_bnd, n_neighs, sys, dx, dt, gm, riemann_solver_flag)
!==============================================================================================================================================|

   implicit none

   integer(4), intent(in) :: nx, ng, n_neighs, sys, riemann_solver_flag
   real(8), intent(in) :: dx, dt, gm
   real(8), dimension(n_neighs, sys), intent(inout) :: flux_bnd
   real(8), dimension(nx), intent(inout) :: ro, mx, my, mz, bx, by, bz, en

   integer(4) :: i, im1
   real(8), dimension(nx) :: ro_pre, mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre
   real(8), dimension(nx) :: fro, fmx, fmy, fmz, fbx, fby, fbz, fen

   real(8) :: dtdx

!----------------------------------------------------------------------------------------------------------------------------------------------|
   dtdx = dt / dx

   call weno_l (fro, fmx, fmy, fmz, fbx, fby, fbz, fen, ro, mx, my, mz, bx, by, bz, en, nx, gm, riemann_solver_flag)

   do i = 2, nx
      im1 = i - 1
      ro_pre(i) = ro(i) - dtdx * (fro(i) - fro(im1))
      mx_pre(i) = mx(i) - dtdx * (fmx(i) - fmx(im1))
      my_pre(i) = my(i) - dtdx * (fmy(i) - fmy(im1))
      mz_pre(i) = mz(i) - dtdx * (fmz(i) - fmz(im1))
      bx_pre(i) = bx(i) - dtdx * (fbx(i) - fbx(im1))
      by_pre(i) = by(i) - dtdx * (fby(i) - fby(im1))
      bz_pre(i) = bz(i) - dtdx * (fbz(i) - fbz(im1))
      en_pre(i) = en(i) - dtdx * (fen(i) - fen(im1))
   enddo

   call weno_l (fro, fmx, fmy, fmz, fbx, fby, fbz, fen, ro_pre, mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre, nx, gm,    &
      riemann_solver_flag)

   do i = 2, nx
      im1 = i - 1
      ro(i) = (ro(i) + ro_pre(i) - dtdx * (fro(i) - fro(im1))) * 0.5d0
      mx(i) = (mx(i) + mx_pre(i) - dtdx * (fmx(i) - fmx(im1))) * 0.5d0
      my(i) = (my(i) + my_pre(i) - dtdx * (fmy(i) - fmy(im1))) * 0.5d0
      mz(i) = (mz(i) + mz_pre(i) - dtdx * (fmz(i) - fmz(im1))) * 0.5d0
      bx(i) = (bx(i) + bx_pre(i) - dtdx * (fbx(i) - fbx(im1))) * 0.5d0
      by(i) = (by(i) + by_pre(i) - dtdx * (fby(i) - fby(im1))) * 0.5d0
      bz(i) = (bz(i) + bz_pre(i) - dtdx * (fbz(i) - fbz(im1))) * 0.5d0
      en(i) = (en(i) + en_pre(i) - dtdx * (fen(i) - fen(im1))) * 0.5d0
   enddo

   flux_bnd(1, 1) = fro(ng)
   flux_bnd(2, 1) = fro(nx - ng)
   flux_bnd(1, 2) = fmx(ng)
   flux_bnd(2, 2) = fmx(nx - ng)
   flux_bnd(1, 3) = fmy(ng)
   flux_bnd(2, 3) = fmy(nx - ng)
   flux_bnd(1, 4) = fmz(ng)
   flux_bnd(2, 4) = fmz(nx - ng)
   flux_bnd(1, 5) = fbx(ng)
   flux_bnd(2, 5) = fbx(nx - ng)
   flux_bnd(1, 6) = fby(ng)
   flux_bnd(2, 6) = fby(nx - ng)
   flux_bnd(1, 7) = fbz(ng)
   flux_bnd(2, 7) = fbz(nx - ng)
   flux_bnd(1, 8) = fen(ng)
   flux_bnd(2, 8) = fen(nx - ng)

!----------------------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine weno
